
# Install and Load Libraries

install.packages(c("ggplot2","psych", "dplyr", "readr", "janitor", "stargazer", "lmtest", "sandwich"))

library(ggplot2)
library(psych)
library(dplyr)
library(readr)
library(janitor)
library(stargazer)
library(lmtest)
library(sandwich)



set.seed(89)


#######################
#### Load and rename datasets 
#######################

Data <- read_csv("Raw_Data.csv")  #Make sure that the CSV file is saved in the same directory as your project to enable easy loading. 
Data = subset(Data, Progress == 100 & Consent == 1)




#######################
#### Data Cleaning
#######################

#Compute age
Data$Year_Born <- as.numeric(Data$Year_Born)
Data$age = 2020 - Data$Year_Born



#Recode gender
Data$Male <- ifelse(Data$Gender == 2, 1, 0)
Data$Male = factor(Data$Male, levels = c(0, 1), labels = c('Female', 'Male'))



#Recode Political position
Data$Political_Position = factor(Data$Political_Position, levels = c(1,2,3,4,5,6), labels = c('Very Left', 'Left', 'Slightly Left', 'Slightly Right', 'Right', 'Very Right'))

Data$Political_Right <- ifelse(Data$Political_Position =="Slightly Right" | Data$Political_Position =="Right" | Data$Political_Position =="Very Right", 1, 0)
Data$Political_Right = factor(Data$Political_Right, levels = c(0, 1), labels = c('Left', 'Right'))


### Recoding Computer Literacy ###
Data$DigitalLiteracy1[is.na(Data$DigitalLiteracy1)] <- 999
Data$DigitalLiteracy2[is.na(Data$DigitalLiteracy2)] <- 999
Data$DigitalLiteracy3[is.na(Data$DigitalLiteracy3)] <- 999
Data$DigitalLiteracy4[is.na(Data$DigitalLiteracy4)] <- 999

answers = Data[c(3:6)]
answers$DigitalLiteracy1 = ifelse(answers$DigitalLiteracy1 == 4,1,0)
answers$DigitalLiteracy2 = ifelse(answers$DigitalLiteracy2 == 1,1,0)
answers$DigitalLiteracy3 = ifelse(answers$DigitalLiteracy3 == 4,1,0)
answers$DigitalLiteracy4 = ifelse(answers$DigitalLiteracy4 == 1,1,0)
answers$score = (answers$DigitalLiteracy1 + answers$DigitalLiteracy2 + answers$DigitalLiteracy3 + answers$DigitalLiteracy4)/4
Data$CompLiteracy = answers$score

tabyl(Data$CompLiteracy)

Data$CompLiteracyBinary <- ifelse(Data$CompLiteracy == 1, 1, 0)
Data$CompLiteracyBinary = factor(Data$CompLiteracyBinary, levels = c(0, 1), labels = c('Low Computer Literacy', 'High Computer Literacy'))


# Compute Education - Binary variable - 1 = academic degree, 0 = no academic degree#
Data$EducationBinary <- ifelse(Data$Education > 5.5, 1, 0)
Data$EducationBinary = factor(Data$EducationBinary, levels = c(0, 1), labels = c('No Academic Degree', 'Academic Degree'))



# Compute Income - Binary variable - 1 = academic degree, 0 = no academic degree#
Data$IncomeBinary <- ifelse(Data$Income > 4.5, 1, 0)
Data$IncomeBinary = factor(Data$IncomeBinary, levels = c(0, 1), labels = c('Lower than average', 'Higher than average'))



# Compute Military Service - Binary variable - 1 = Served, 0 = Hasn't served #
Data$MilitaryService[Data$MilitaryService == 2] <- 0
Data$MilitaryService = factor(Data$MilitaryService, levels = c(0, 1), labels = c('Didnt Serve', 'Served'))


# Compute Employment - Binary variable - 1 = Employed, 0 = Unemployed #
Data$EmploymentBinary <- ifelse(Data$Employment < 3.5, 1, 0)
Data$EmploymentBinary = factor(Data$EmploymentBinary, levels = c(0, 1), labels = c('Not Employed', 'Employed'))




#### Military assertiveness
Data$Military_Assert_1 = as.numeric(Data$Military_Assert_1)
Data$Military_Assert_2 = as.numeric(Data$Military_Assert_2)
Data$Military_Assert_3 = as.numeric(Data$Military_Assert_3)

Data$Military_Assert_2_re[Data$Military_Assert_2 == 1] <- 5
Data$Military_Assert_2_re[Data$Military_Assert_2 == 2] <- 4
Data$Military_Assert_2_re[Data$Military_Assert_2 == 3] <- 3
Data$Military_Assert_2_re[Data$Military_Assert_2 == 4] <- 2
Data$Military_Assert_2_re[Data$Military_Assert_2 == 5] <- 1

Data$MilAssert <- (Data$Military_Assert_1 + Data$Military_Assert_2_re + Data$Military_Assert_3)/3

Data$MilAssertBinary <- ifelse(Data$MilAssert > 3.1, 1, 0)
Data$MilAssertBinary = factor(Data$MilAssertBinary, levels = c(0, 1), labels = c('Low Military Assertiveness', 'High Military Assertiveness'))



#International trust
Data$Int_Trust = factor(Data$Int_Trust, levels = c(1, 2), labels = 
                          c('High Trust', 'Low Trust'))




#National chauvinism
Data$Nat_Chauv1 = as.numeric(Data$Nat_Chauv1)
Data$Nat_Chauv2 = as.numeric(Data$Nat_Chauv2)

Data$Nat_Chauv1_re[Data$Nat_Chauv1 == 1] <- 4
Data$Nat_Chauv1_re[Data$Nat_Chauv1 == 2] <- 3
Data$Nat_Chauv1_re[Data$Nat_Chauv1 == 3] <- 2
Data$Nat_Chauv1_re[Data$Nat_Chauv1 == 4] <- 1

Data$NatChauv <- (Data$Nat_Chauv1_re + Data$Nat_Chauv2)/2

tabyl(Data$NatChauv)



Data$AttentionCheck1
Data$ConditionCode

### Computing Attention Check ###
Data$attention1[Data$ConditionCode == 1 & Data$AttentionCheck1 == 1] <- 1
Data$attention1[Data$ConditionCode == 2 & Data$AttentionCheck1 == 1] <- 1
Data$attention1[Data$ConditionCode == 3 & Data$AttentionCheck1 == 1] <- 1
Data$attention1[Data$ConditionCode == 4 & Data$AttentionCheck1 == 1] <- 1
Data$attention1[Data$ConditionCode == 5 & Data$AttentionCheck1 == 1] <- 1

Data$attention2[Data$ConditionCode == 1 & Data$AttentionCheck2 == 2] <- 1
Data$attention2[Data$ConditionCode == 2 & Data$AttentionCheck2 == 1] <- 1
Data$attention2[Data$ConditionCode == 3 & Data$AttentionCheck2 == 1] <- 1
Data$attention2[Data$ConditionCode == 4 & Data$AttentionCheck2 == 1] <- 1
Data$attention2[Data$ConditionCode == 5 & Data$AttentionCheck2 == 1] <- 1

Data$attention3[Data$ConditionCode == 1 & Data$AttentionCheck3 == 1] <- 1
Data$attention3[Data$ConditionCode == 2 & Data$AttentionCheck3 == 1] <- 1
Data$attention3[Data$ConditionCode == 3 & Data$AttentionCheck3 == 2] <- 1
Data$attention3[Data$ConditionCode == 4 & Data$AttentionCheck3 == 3] <- 1
Data$attention3[Data$ConditionCode == 5 & Data$AttentionCheck3 == 3] <- 1

Data$attention1[is.na(Data$attention1)] <- 0
Data$attention2[is.na(Data$attention2)] <- 0
Data$attention3[is.na(Data$attention3)] <- 0

Data$AttentionCheckTotal = (Data$attention1 + Data$attention2 + Data$attention3)

tabyl(Data$AttentionCheckTotal)




### Recoding main dependent variable into 7 point scale ###
Data$Approval[Data$Pres_Approval_App == 1] <- 7
Data$Approval[Data$Pres_Approval_App == 2] <- 6
Data$Approval[Data$Pres_Approval_Dis == 1] <- 1
Data$Approval[Data$Pres_Approval_Dis == 2] <- 2
Data$Approval[Data$Pres_Approval_Lean == 1] <- 5
Data$Approval[Data$Pres_Approval_Lean == 2] <- 3
Data$Approval[Data$Pres_Approval_Lean == 3] <- 4

Data$Approval = as.numeric(Data$Approval, levels = c(1, 2, 3, 4, 5, 6, 7), labels = 
                             c('Strongly Disapprove', 'Disapprove', 'Slightly Disapprove', 'Neutral', 'Slightly Approve', 'Approve', 'Strongly Approve'))





### Recoding secondary dependent variable (approval of Biden performance) into 7 point scale ###
Data$BidenApproval[Data$Situational_Appr_2 == 1] <- 7
Data$BidenApproval[Data$Situational_Appr_2 == 2] <- 6
Data$BidenApproval[Data$Situational_Appr_3 == 1] <- 1
Data$BidenApproval[Data$Situational_Appr_3 == 2] <- 2
Data$BidenApproval[Data$Situational_Appr_4 == 1] <- 5
Data$BidenApproval[Data$Situational_Appr_4 == 2] <- 3
Data$BidenApproval[Data$Situational_Appr_4 == 3] <- 4

Data$BidenApproval = as.numeric(Data$BidenApproval, levels = c(1, 2, 3, 4, 5, 6, 7), labels = 
                                  c('Strongly Disapprove', 'Disapprove', 'Slightly Disapprove', 'Neutral', 'Slightly Approve', 'Approve', 'Strongly Approve'))








### Generate dummy variables for treatment conditions ###
      #1: Stay Out, 2: Engage, 3: Not Engage, 4: Cyber Degradation, 5: Cyber Disruption

      #This variable includes all five conditions
Data$Strategy = factor(Data$ConditionCode, levels = c(1, 2, 3, 4, 5), labels = 
                         c('Stay Out', 'Back Down', 'Engage', 'Engage Cyber Major', 'Engage Cyber Minor'))

      #This variable combines the two cyber treatments into a single group
Data$Strategy4 = factor(Data$ConditionCode)
Data$Strategy4[Data$ConditionCode == 1] <- 1
Data$Strategy4[Data$ConditionCode == 2] <- 2
Data$Strategy4[Data$ConditionCode == 3] <- 3
Data$Strategy4[Data$ConditionCode == 4] <- 4
Data$Strategy4[Data$ConditionCode == 5] <- 4

Data$Strategy4 = factor(Data$Strategy4, levels = c(1, 2, 3, 4), labels = 
                          c('Stay Out', 'Back Down', 'Engage Militarily', 'Engage Cyber'))



    #Dummy Variables
Data$StayOut <- ifelse(Data$Strategy=="Stay Out",1,0)
Data$BackDown <- ifelse(Data$Strategy=="Back Down",1,0)
Data$Engage <- ifelse(Data$Strategy=="Engage",1,0)
Data$CyberDegradation <- ifelse(Data$Strategy=="Engage Cyber Major",1,0)
Data$CyberDisruption <- ifelse(Data$Strategy=="Engage Cyber Minor",1,0)
Data$CyberAll <- ifelse(Data$Strategy=="Engage Cyber Minor" | Data$Strategy=="Engage Cyber Major",1,0)




##################################
########    Figure 2    ########  
#################################

m1a_approval<- lm(Approval ~ BackDown + Engage + CyberAll, data=Data) 
summary(m1a_approval)
out_approval <- coeftest(m1a_approval, vcov = vcovHC(x = m1a_approval, type = "HC2")) # For robust standard errors 
ggdf2 <- data.frame(treatment = c("Stays out", "Backs Down", "Engage Militarily", "Cyber Attack"),
                    est = c(0, out_approval[2:4, 1]),
                    se = c(NA, out_approval[2:4, 2]))
ggdf2$treatment <- factor(ggdf2$treatment, levels = unique(ggdf2$treatment))

#png("Figure2.png", res=600, height=10, width=10, units="in")     # Add this line and the dev.off() line to save the figure
ggplot(ggdf2, aes(x = treatment, y = est)) + geom_point() +
  geom_errorbar(aes(ymin = est - 1.96 *se, ymax = est + 1.96 * se), width = .2) +
  theme_minimal() + xlab("\nExperimental condition") + ylab("Average Treatment Effect\n") +
  theme(                                                     # this block removes the vertical gridlines
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank(),
    text = element_text(size=16)
  ) +
  geom_hline(yintercept = 0, lty = 3, col = "red")
  # + scale_x_discrete(guide = guide_axis(n.dodge=2))         # This would puts the x-axis labels on 2 lines to avoid overlapping
#dev.off()





##################################
########    Figure 3    ########  
#################################

m2_approval_hte <- lm(Approval ~ CompLiteracyBinary*(BackDown + Engage + CyberDegradation + CyberDisruption), 
                      data=Data) 
out_hte2 <- coeftest(m2_approval_hte, vcov = vcovHC(x = m2_approval_hte, type = "HC2")) # For robust standard errors 
vcov_hte2 <- vcovHC(x = m2_approval_hte, type = "HC2")
se_3 <- function(row1, row2, row3, vcov){
  return(sqrt(vcov[row1, row1] + vcov[row2, row2] + vcov[row3, row3] + 
                2 * vcov[row1, row2] + 2 * vcov[row2, row3] + 2 * vcov[row1, row3]))
}
rownames(vcov_hte2)
ggdf3 <- data.frame(treatment = rep(c("Stays out", "Backs Down", "Engage Militarily", 
                                      "Degradative Cyber Attack", "Disruptive Cyber Attack"), 2),
                    Type = rep(c("Low Tech Literacy", "High Tech Literacy"), each = 5),
                    est = c(0, out_hte2[3:6, 1], out_hte2[2, 1], 
                            out_hte2[2, 1] + out_hte2[3:6, 1] + out_hte2[7:10, 1]),
                    se = c(NA, out_hte2[3:6, 2], out_hte2[2, 2], 
                           se_3(2, 3, 7, vcov = vcov_hte2), 
                           se_3(2, 4, 8, vcov = vcov_hte2),
                           se_3(2, 5, 9, vcov = vcov_hte2),
                           se_3(2, 6, 10, vcov = vcov_hte2)))
ggdf3$treatment <- factor(ggdf3$treatment, levels = unique(ggdf3$treatment))
ggdf3$Type <- factor(ggdf3$Type, levels = unique(ggdf3$Type))


#png("Figure3", res=600, height=10, width=10, units="in")         # Add this line and the dev.off() line to save the figure
ggplot(ggdf3, aes(x = treatment, y = est, col = Type)) + 
  geom_point(aes(shape=Type), position = position_dodge(width = 0.2)) +
  geom_errorbar(aes(ymin = est - 1.96 *se, ymax = est + 1.96 * se), width = .2,
                position = position_dodge(width = 0.2)) +
  theme_minimal() + xlab("\nExperimental Condition") + ylab("Average Treatment Effect\n") + 
  theme(                                           # these lines in the bracket removes the vertical gridlines and changes font size
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank(),
    text = element_text(size=16),
    axis.text.x  = element_text(angle=45, vjust=0.5)         # This turns the axis text on an angle
  ) +
  geom_hline(yintercept = 0, lty = 3, col = "red") 
#dev.off()





##################################
########    Figure 4    ########  
#################################

m1_approval_mil_assert <- lm(Approval ~ MilAssertBinary *(BackDown + Engage + CyberDegradation + CyberDisruption), 
                             data=Data) 
out_hte <- coeftest(m1_approval_mil_assert, vcov = vcovHC(x = m1_approval_mil_assert, type = "HC2")) # For robust standard errors 
vcov_hte <- vcovHC(x = m1_approval_mil_assert, type = "HC2")
se_3 <- function(row1, row2, row3, vcov){
  return(sqrt(vcov[row1, row1] + vcov[row2, row2] + vcov[row3, row3] + 
                2 * vcov[row1, row2] + 2 * vcov[row2, row3] + 2 * vcov[row1, row3]))
}
rownames(vcov_hte)
ggdf2 <- data.frame(treatment = rep(c("Stays out", "Backs Down", "Engage Militarily", 
                                      "Degradative Cyber Attack", "Disruptive Cyber Attack"), 2),
                    Type = rep(c("Doves", "Hawks"), each = 5),
                    est = c(0, out_hte[3:6, 1], out_hte[2, 1], 
                            out_hte[2, 1] + out_hte[3:6, 1] + out_hte[7:10, 1]),
                    se = c(NA, out_hte[3:6, 2], out_hte[2, 2], 
                           se_3(2, 3, 7, vcov = vcov_hte), 
                           se_3(2, 4, 8, vcov = vcov_hte),
                           se_3(2, 5, 9, vcov = vcov_hte),
                           se_3(2, 6, 10, vcov = vcov_hte)))
ggdf2$treatment <- factor(ggdf2$treatment, levels = unique(ggdf2$treatment))
ggdf2$Type <- factor(ggdf2$Type, levels = unique(ggdf2$Type))


#png("Figure4", res=600, height=10, width=10, units="in")
ggplot(ggdf2, aes(x = treatment, y = est, col = Type)) + 
  geom_point(aes(shape=Type), position = position_dodge(width = 0.2)) +
  geom_errorbar(aes(ymin = est - 1.96 *se, ymax = est + 1.96 * se), width = .2,
                position = position_dodge(width = 0.2)) +
  theme_minimal() + xlab("\nExperimental Condition") + ylab("Average Treatment Effect\n") +
  theme(                                           # these lines in the bracket removes the vertical gridlines and changes font size
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank(),
    text = element_text(size=16),
    axis.text.x  = element_text(angle=45, vjust=0.5)         # This turns the axis text on an angle
  ) +
  geom_hline(yintercept = 0, lty = 3, col = "red") 
# + scale_x_discrete(guide = guide_axis(n.dodge=2))         # This puts the labels on 2 lines to avoid overlapping
#dev.off()





##################################
########    Appendix B - Table A1    ########  
#################################

get_mean_sd <- function(vector){
  mean <- mean(vector, na.rm = T)
  sd <- sd(vector, na.rm = T)
  sem <- sd/(sqrt(length(vector)))
  output <- paste(round(mean, 3), " (", round(sem, 3), ")", collapse = "", sep = "")
  return(output)
}

balance_data <- group_by(Data, Strategy) %>%
  summarize(Age = get_mean_sd(age),
            Gender = get_mean_sd(as.numeric(Male)), 
            PoliticalIdeology = get_mean_sd(as.numeric(Political_Position)),
            EducationBinary = get_mean_sd(as.numeric(EducationBinary)),
            Income = get_mean_sd(as.numeric(Income)),
            MilitaryAssertiveness = get_mean_sd(as.numeric(MilAssert)))
balance_data

#The following code will export the balance table to a word document
#out <- t(balance_data[,2:7])
#colnames(out) <- c("Stays Out", "Backs Down", "Engages", "Engages Cyber Degradation", "Engages Cyber Disruption")
#stargazer(out)
#out



##################################
########    Appendix B - Table A2    ########  
#################################

ProbitBalanceCheck_StayOut <- glm(StayOut ~  age + Male + EducationBinary + IncomeBinary + Political_Right + 
                                    MilAssertBinary + EmploymentBinary + MilitaryService, data=Data) 
summary(ProbitBalanceCheck_StayOut)
logLik(ProbitBalanceCheck_StayOut)
AIC(ProbitBalanceCheck_StayOut)
ProbitBalanceCheck_BackDown <- glm(BackDown ~  age + Male + EducationBinary + IncomeBinary + Political_Right +
                                     MilAssertBinary + EmploymentBinary + MilitaryService, data=Data) 
summary(ProbitBalanceCheck_BackDown)
logLik(ProbitBalanceCheck_BackDown)
AIC(ProbitBalanceCheck_BackDown)
ProbitBalanceCheck_Engage <- glm(Engage ~  age + Male + EducationBinary + IncomeBinary + Political_Right +
                                   MilAssertBinary + EmploymentBinary + MilitaryService, data=Data) 
summary(ProbitBalanceCheck_Engage)
logLik(ProbitBalanceCheck_Engage)
AIC(ProbitBalanceCheck_Engage)
ProbitBalanceCheck_CyberDegradation <- glm(CyberDegradation ~  age + Male + EducationBinary + IncomeBinary + Political_Right +
                                             MilAssertBinary + EmploymentBinary + MilitaryService, data=Data) 
summary(ProbitBalanceCheck_CyberDegradation)
logLik(ProbitBalanceCheck_CyberDegradation)
AIC(ProbitBalanceCheck_CyberDegradation)
ProbitBalanceCheck_CyberDisruption <- glm(CyberDisruption ~  age + Male + EducationBinary + IncomeBinary + Political_Right +
                                            MilAssertBinary + EmploymentBinary + MilitaryService, data=Data) 
summary(ProbitBalanceCheck_CyberDisruption)
logLik(ProbitBalanceCheck_CyberDisruption)
AIC(ProbitBalanceCheck_CyberDisruption)

#The code below exports the five regression plots above into a word document
#         stargazer(ProbitBalanceCheck_StayOut, ProbitBalanceCheck_BackDown, ProbitBalanceCheck_Engage, 
#          ProbitBalanceCheck_CyberDegradation, ProbitBalanceCheck_CyberDisruption,
#          title="Probit Regression", align=TRUE,
#          type = "html",
#          star.cutoffs = c(0.05, 0.01, 0.001),       # Add this line into stargazer output to change the significance values that are associated with stars 
#          out = "star_descriptive_2.doc")





##################################
########    Appendix C - Table A3    ########  
#################################

columnA <- lm(BidenApproval ~ BackDown + Engage + CyberAll, data=Data) 
summary(columnA)
columnB <- lm(BidenApproval ~ BackDown + Engage + CyberAll + 
                age + Male + EducationBinary + IncomeBinary + Political_Right +
                CompLiteracyBinary + MilAssertBinary, data=Data) 
summary(columnB)
columnC <- lm(BidenApproval ~ BackDown + Engage + CyberDegradation + CyberDisruption +
                CompLiteracyBinary:BackDown + CompLiteracyBinary:Engage + 
                CompLiteracyBinary:CyberDegradation + CompLiteracyBinary:CyberDisruption,
              data=Data) 
summary(columnC)
columnD <- lm(BidenApproval ~ BackDown + Engage + CyberDegradation + CyberDisruption +
                MilAssertBinary:BackDown + MilAssertBinary:Engage + 
                MilAssertBinary:CyberDegradation + MilAssertBinary:CyberDisruption,
              data=Data) 
summary(columnD)

#The code below exports the five regression plots above into a word document
#           stargazer(columnA, columnB, columnC, columnD,
#          title="columnsAlternative", align=TRUE,
#          type = "html",
#          star.cutoffs = c(0.05, 0.01, 0.001),       # Add this line into stargazer output to change the significance values that are associated with stars 
#          out = "columnsAlternative.doc")

#     star.cutoffs = c(0.05, 0.01, 0.001)     




##################################
########    Appendix D - Attention Checks    ########  
#################################

tabyl(Data, AttentionCheckTotal, Strategy)  # This shows the number of people who answered 0,1, 2 or 3 attention check questions correctly in each treatment group



##################################
########    Appendix E - Table A4    ########  
#################################

column1 <- lm(Approval ~ BackDown + Engage + CyberAll, data=Data) 
summary(column1)
column2 <- lm(Approval ~ BackDown + Engage + CyberAll + 
                age + Male + EducationBinary + IncomeBinary + Political_Right +
                CompLiteracyBinary + MilAssertBinary, data=Data) 
summary(column2)
column3 <- lm(Approval ~ BackDown + Engage + CyberDegradation + CyberDisruption +
                CompLiteracyBinary:BackDown + CompLiteracyBinary:Engage + 
                CompLiteracyBinary:CyberDegradation + CompLiteracyBinary:CyberDisruption,
              data=Data) 
summary(column3)
column4 <- lm(Approval ~ BackDown + Engage + CyberDegradation + CyberDisruption +
                MilAssertBinary:BackDown + MilAssertBinary:Engage + 
                MilAssertBinary:CyberDegradation + MilAssertBinary:CyberDisruption,
              data=Data) 
summary(column4)


#The code below exports the five regression plots above into a word document
#           stargazer(column1, column2, column3, column4,
#          title="columns", align=TRUE,
#          type = "html",
#          star.cutoffs = c(0.05, 0.01, 0.001),       # Add this line into stargazer output to change the significance values that are associated with stars 
#          out = "columns.doc")

